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<**] . Abstract. Asynchronous dynamics given by the master equation in the Sherrington- 

\ Kirkpatrick (SK) spin-glass model is studied based on dynamical replica theory 

(DRT) with an extension to take into account the autocorrelation function. The 
dynamical behaviour of the system is approximately described by dynamical equations 
of the macroscopic quantities: magnetization, energy contributed by randomness, and 
the autocorrelation function. The dynamical equations under the replica symmetry 
assumption are derived by introducing the subshell equipartitioning assumption and 
£h ■ exploiting the replica method. The obtained dynamical equations are compared with 

Monte Carlo (MC) simulations, and it is demonstrated that the proposed formula 
describes well the time evolution of the autocorrelation function in some parameter 
regions. The study offers a reasonable description of the autocorrelation function in 
the SK spin-glass system. 



PACS numbers: 75.10.Nr, 75.40. Gb 
1. Introduction 

Random systems that possess frustration have many metastable states in their free 
energy landscape, which are expected to contribute to the system behaviour. Spin-glass 
systems are representative of such systems and their equilibrium properties have been 
extensively studied in a mean-field model, the Sherrington-Kirkpatrick (SK) model. 
Mean-field theory based on the replica method reveals the existence of a transition 
associated with replica symmetry breaking (RSB) that is interpreted with the aid of 
Thouless-Anderson-Palmer (TAP) theory as the exponential appearance of metastable 
states in the free energy landscape [HE]- A quantitative understanding of the metastable 
states is attempted by counting the number of solutions of the TAP equation or belief 
propagation equation [3j HJ |5] . 

Below a critical temperature, spin glasses in experiments show slow relaxation 
dynamics and an aging effect in which the dynamical properties depend on the history 
even after a long time period [6]. The dynamical behaviour and equilibrium properties 
are considered to be related, but the connection is not fully understood. It has been 
suggested that dynamical behaviour of the SK model is qualitatively similar to spin 
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glasses in experiments even though the realistic spin glasses have finite dimensions 
[7J. Therefore, clarifying the relaxation dynamics in the SK model provides insight 
into the connection between the equilibrium and dynamical properties, as well as into 
the dynamics of random systems. For this purpose, an analytical description of the 
relaxation dynamics in the SK model that can be compared with our knowledge of the 
equilibrium states would be of significant use. 

Dynamical replica theory (DRT) provides a statistical-mechanical description of the 
dynamics in random systems [8]. In DRT, the time evolution of macroscopic quantities 
is derived according to the replica method. A characteristic feature of DRT is the 
introduction of an assumption called subshell equipartitioning in which the microscopic 
states in the subshell that are characterized by a certain macroscopic quantity value 
appear with the same probability. By introducing this assumption, one can obtain closed 
dynamical equations that describe the time evolution of the macroscopic quantities. The 
derivation of the dynamical equations and their DRT solutions are tractable even if the 
relaxation dynamics over long time are considered, although an exact description of the 
dynamics is sacrificed. DRT offers a contrasting perspective compared to the generating 
functional method [9j [TO] in which the dynamical equations are exact but difficult to 
solve over long time steps. 

The physical implications of subshell equipartitioning and its influence on the 
accuracy of the theoretical prediction have been studied [HJ [12] • The assumption 
removes the microscopic time correlation, and the system only depends on the past value 
of the macroscopic quantities that characterize the subshell. It is well known that the 
time correlation is significant in describing the dynamics of random systems. In fact, 
it has been demonstrated that analytical descriptions of the dynamics are improved 
by taking into account the time correlation in the Hopfield model [131 EH] an d in the 
interference canceller of code-division-multiple-access (CDMA) multiuser detection [15]. 
However, DRT in its present formalism cannot deal with the time correlation, and hence 
to apply DRT to random systems and discuss the time correlation, an extension should 
be introduced. 

A central quantity in the study of the dynamics of spin glasses is the spin 
autocorrelation function that describes the correlation between microscopic states at two 
time steps. It is experimentally and numerically observed with the conjugate response 
function as a quantity that indicates the slow relaxation dynamics [6] and its analytical 
description has been studied in particular Langevin dynamics of soft spins [HI [17] . In 
this paper, I study the sequential dynamics of the SK model with DRT while including 
the time evolution of the autocorrelation function. DRT has been studied mainly for two 
cases where the subshell is characterized by a magnetization and energy set contributed 
by randomness [8] or a conditional local field distribution [18]. I start with the simplest 
DRT [8] and introduce the autocorrelation function into the formulation. The evolution 
equation of the autocorrelation function is derived based on the replica method and its 
dependence on the phase is discussed. 

This paper is organized as follows. In section [2j I explain the model settings 
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and time evolution of the joint probability of microscopic states at two time steps, 
which is described by Glauber dynamics [15] . In section [3j macroscopic quantities, 
magnetization, energy contributed by randomness, and the autocorrelation function, 
are introduced to characterize the microscopic state space. The closed formula of the 
evolutional dynamics of the quantities is obtained following the procedures of DRT, and 
the expressions under the replica symmetry (RS) assumption are derived. In section HI 
the obtained dynamical equations are numerically solved and the results are compared 
with MC simulations. Finally, section [5] is devoted to conclusions and the outlook for 
further developments. 

2. Model setup 

The system being studied here is a mean-field spin-glass model, the SK model, which 
consists of N Ising spins S = {Si} (i = 1, ■ ■ ■ , N). The Hamiltonian is given by 

H(S\J) = -£ JySiSj -h^ S u (1) 

i<3 i 

where J^- is the interaction between spins Si and Sj, and h is the external field. The 
interactions J = {Jij} are symmetric (J^ = Jji) and assumed to be independently and 
identically distributed according to the Gauss distribution with a mean of Jq/N and a 
variance of J 2 /N; 

W,0 = y / ^exp[-A( Jy ._^) 2 ]. (W) (2) 

Here, I consider the conditional probability of a microscopic spin configuration 
cr at t w + t under a given spin configuration s at time t w and interaction J as 
p(er;t w + t\s;t w ,J). The time t w is called the waiting time. The two configurations 
cr and s are related to each other through cr = s o d, where d = {oli} (i = 1, • • • , N) 
is a vector consisting of Ising variables, and o denotes the Hadamard product, which is 
a product with respect to each component; the i-th component of s o d is Sjdj. When 
di = 1, the i-th spin configuration at time t w and that at time t w + t are the same; 
otherwise they differ. The time evolution of the microscopic conditional probability in 
the direction of t is given by Glauber dynamics, which is an asynchronous update of 
spin configurations described by the master equation, 

d N 

— p(s od;t w + t\s; t w , J) = ^{p{s o F k d; t w + t\s; t w , J)w k (s o F k d\J) 

k=l 

- p(s od;t w + t\s; t w , J)w k (s o d\J)}, (3) 

where F k is an operator that flips the configuration of the k-th Ising variable as 
F k s = {si, ■ ■ ■ , — s k , ■ ■ ■ , Sjv}. The transition probability under a given J, w k (s\J), 
is given by 

Wk(s\J) = -(1 - s fe tanh(/3/i fc (s|J))), (4) 
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where h k {s\J) = Y^i^k JikSi + h is the local field of the k-th spin. The time evolution 
of the joint probability p(cr; t w + r, s; t w \J) is obtained by multiplying both sides of (j3J) 
by the probability distribution of a microscopic state at time t w , p(s;t w \J). The fixed 
point of ([3]) that is obtained at t — > oo corresponds to the equilibrium distribution 

p(cr; t w + t, s; t w \J) ->■ exp(-/3F(o-| J)), (5) 

where Z is the normalization constant. 



3. Time evolution of macroscopic quantities 

The macroscopic states are defined by the following macroscopic quantities, 

1 N 1 J 

m o(s) = r o( s ) = T7 E ( J *i - 7^) s ^i ( 6 ) 

j AT ^ j 

■y N ^ N 

iv i=l iv i=l 

where mo and ro are the magnetization and the absolute value of the energy contributed 
by randomness (referred to as randomness energy, hereafter) at time t w , respectively, and 
m and r are the corresponding values at time t w +t. The quantity c is the autocorrelation 
function between time t w +t and t w , and it is a function of d. The probability distribution 
of the macroscopic quantities is given by 

Pt w +t,t w (^\ J ) = E^( cr ' t «' +t,s;t w \J)5(m -m(s))5(r -r(s)) 

<T,S 

x 8{m - m{cr))8{r - r(cr))8{c - c(d)), (9) 

where f2 = {mo, ro, m, r, c}. There are three Ising variables cr, s, and d, but the 
relationship cr = s o d reduces the number of the independent variables to two. The 
double summation over the microscopic states of a and s should be taken in mind of 
the relationship cr = s o d. The time evolution of P tw+ t ) t w {^\ J) in the direction of t is 
derived by substituting (j3J) and @ into (jUJ) and employing the Taylor expansion with 
respect to 1/N: 

-^Pt w +t,t w (Q\J) 

d 1 N 



a 1 

- ^-P tlu+tli Jfi| J)(- E((tanh/3(J m + J* fc (<r| J) + h))) n - m) 

fc=i 

d I N 

- J)(^ ^((^tanh/3(J m + Jz k (cr\J) + h)}) n - 2< 

k=l 

d I N 

- Q-Pt w +t,t w (n\J)(jj £ ((s fc tanh/3( J m + Jz k (cr\J) + h))) a - cpO) 

k=l 
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where z k (er\ J) = Y^i^ki^ki ~ Jo/N)^ is the local field contributed by randomness of the 
fc-th spin. The notation ((■ ■ -}}^ represents the average within the subshell specified by 
the set of macroscopic quantities Q as 

= S«r, a • • t w + t, s; t w \J)W(tr, s\Q) 
{{ "' ))n E<r, a p{<r;t w + t,s;t w \J)W{<T,s\n) ' 1 ' 

where W(cr, s\Q) = S(m — m (s))5(r — r (s))8(m — m(cr))5(r — r(cr))S(c — c(d)). By 
regarding (TlOT) as the total differential form with respect to m, r, and c, their dynamical 
equations are given by 
d r 

—m= dz D{z,S\Q, J) tanh(/3( J m + Jz + h)) - m (12) 



S=±l 



d 
eft 
d 

— c 
dt 



r = f dz D(z,S\n, J) tanh(/3( J m + Jz + h)) -It (13) 
^ 5=±l 

/ dz ^ D^S'in, J)5 tanh(/3(J m+ Jz + h)) - c, (14) 
^ s=±i 

where S takes the value ±1, and z is the effective noise caused by quenched randomness. 
The effective noise is distributed according to D(z, S\Q, J), which is given by 

D(z, S\Q, J) = lim (15 ) 

where 5 Skt s is the Kronecker delta. 

3.1. Closed formula of the dynamical equations 

Following the procedures of DRT, a closed formula of the dynamical equations of the 
macroscopic variables is obtained by assuming the realization of two properties: self- 
averaging and subshell equipartitioning. The self-averaging property guarantees that 
the noise distribution, D(z, S\Q, J), which depends on a realization of J, converges to a 
typical distribution D(z, S\fl) = [D(z, S\fl, J)]j in the limit N — > <x>, where [• • ■] j is the 
average over the randomness J according to the probability distribution (j2j). Subshell 
equipartitioning ensures that the microscopic states within a subshell specified by the 
macroscopic quantities Q = {m Q ,r Q ,m,r,c} appear with the same probability at each 
time step. As a consequence of these assumptions, the noise distribution D(z, S\fl, J) 
is replaced with D(z, S\Q) as given by 

D{z,bW-hmJ Ea>s W(<r, S \n) \f (16) 

At equilibrium, the probability distribution of the microscopic states is given by 
the Boltzmann distribution with a Hamiltonian ([T|) that can be expressed in terms of 
m and r, and hence the subshell equipartitioning property with m and r is realized. 
However, the property is not expected to hold in cases far from equilibrium, and to 
apply the subshell equipartitioning assumption to non-equilibrium states is crucial for 
the theoretical accuracy. 
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The average over the quenched randomness in fll6p is calculated by introducing the 
replica expression as 

1 N 

D(,,^) = hmjim £ £ [^E^-^(^))^ 

a x ,—,tT n a 1 ,— ,s n fc=l 
n 

x n^slfil- (17) 

a=l J 

The right-hand side of (1TT1) is calculated for a positive integer n and then analytically 
continued to a non-integer n by taking the limit to 0. After the calculations, the noise 
distribution can be expressed as 

D(z,S\Q) = lim lim -Le^- p (^ 5 exp{-^ - z) 2 }h 
where P = \ E ab (p a p b + PoP b + p b p a + P0P0), z = pV + p a s a c + E£ =2 (pVY a + Po s a q} a ), 

$ = extr [ - £(m/i a + m /ig + rp a + r p a + c 7 a ) + - £(p a2 + 2c 2 pVS + pf 



a=l a 



m 



~ ~ E(P W + PolfPo + Po<?f V + rfrf 8 V) + log £ e S 
S = £(pV a + P> a + 7V a s a ) 

a 

+ £(<f V a pW + g^PoPo 5 ^ 6 + q?p a oP b s a ° b + q ba p b p a s b (r a ), (19) 

and (A) s = E{^}>{^> ^ eH / £{<r«},{-»} e H The variables ft = {{pg}, {pg}, {p a }, {p a }, {7°}} 
and Q = {{qo b }, {q ab }, {q? b }} (a, & = 1, •••,n, a 7^ b) are the conjugates of fl — 
{mo, r , m, r, c} and the overlap parameters, respectively. The notation extr^ Q .F(f2, Q) 
denotes the extremization of a function Q) with respect to Q and Q. At the 

extremum of ( TTBj) . the variables are related to the physical quantities through 

-E^t —1 m = -}_^—^ g-, (20) 

i- £(g ab V + <?f 2 Po), r„ = i E(% afe2 Po + sfV) (21) 

aft ab 

1 ^ E CT s a a s a e s , . 

-E y e s . ( 22 ) 

E<T,S ^ E(T,S ^ E<T,S ^ 

3.2. Replica symmetry 

I introduce the RS assumption here for the conjugate and overlap parameters as 

p a = p, pg = Po, 7 a = 7, (24) 
q ab = q, qo" = Qo, qf = ft, (25) 
p a = p, pg = p . (26) 



and 
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With these assumptions, (I2"0"j) -( l2"3"j) are transformed to 

tanh X + tanh X tanh Y 



m 



DuDv 



DuDv 



DuDv 



1 + tanh Xo tanh X tanh Y 
tanh X + tanh X tanh Y 

1 + tanh Xo tanh X tanh Y 
tanh Y + tanh Xo tanh X 
1 + tanh Xq tanh X tanh Y 



P 



Po 



2{r (c 2 -g 2 )-r(l-g 2 )} 

(c 2 - g 2 ) 2 - (1 - g 2 )(l - Q 2 )' 
2{r(c 2 -g 2 )-r (l-g 2 )} 

(c 2 -g 2 ) 2 -(l-g 2 )(l-g 2 )' 



and 



DuDv(^- 



tanh X + tanh Xq tanh Y 



qo 



q t 



DuDv 



DuDv 



+ tanh X tanh X tanh Y 
tanh X + tanh X tanh K 
1 + tanh X tanh X tanh Y 
tanh X + tanh X tanh Y 

1 + tanh X tanh X tanh Y 1 + tanh X tanh X tanh K 



tanh Xq + tanh X tanh Y 



(27) 
(28) 
(29) 

(30) 
(31) 

(32) 
(33) 
(34) 



where 



X = p + pu^fq 



X = fi + p M' 



\ 



go 



g 



(35) 

(36) 
(37) 



y = 7 - gtppo- 

The noise distribution under the RS assumption is given by a function of the RS 
quantities as 



D(z,S\fl) 



DuDv 



8V27T 



x{(l + S) 

+ 0--S)- 



[1 + tanh Xi)(l + tanh X 1 °)(l + tanh Y) 



1 + tanh Xi tanh X% tanh K 



tanhX 2 )(l -tanhX 2 °)(l 



tanhy) _i 
— e 2 



+ (1 + 5)- 
+ (1-5)- 



1 + tanh X 2 tanh X^ tanh Y 
[1 - tanhX 3 )(l + tanhX 3 °)(l - tanhY) 
1 + tanh X 3 tanh X$ tanh Y ' 
;i -tanhX 4 )(l -tanhX 4 °)(l + tanhF) 



(z-z 2 f 



i(z-Z4,) 



(38) 



1 + tanh X 4 tanh X 4 tanh Y 
which is a mixture of four Gaussian distributions. The means of the four Gaussian 
distributions are given by z\ = p(l — q) + po(c — q t ), ~z 2 = p(l — q) — Po(c — qt), 
z 3 = — p(l — g) + p (c — qt), and 24 = — p(l — g) — po(c — g t ), and the functions Xj and 
X 4 ° for z = 1, • ■ ■ , 4 are given by 



Xi = n + py/q{ \Jl - qu + y/q(z - z { )}, 



(39) 
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{ vl {y/ 1 ~ qu + ^ z ~ Zi> >) + \ 




(40) 



3. 3. Fixed point of the macroscopic equations 

The dynamical equations reach the fixed point, for which dm/dt = 0, dr/dt = 0, and 
dc/dt = 0, after a long time evolution of t — > oo. At the fixed point denoted by m*, r*, 
and c*, the following relationships hold: 

m* = I Dz lim (tanh (3(J m + J(z + z) + h))* (41) 
r* =- J Dz \im ((z + z) tanh (3(J m + J(z + z) + h)) E (42) 

c* = f Dz lim (s 1 tanh /3(J m + J(z + z) + h)) s , (43) 

where f Dz = (27T)" 1 / 2 dz exp(-z 2 /2). Equations (jiB- fpjl are obtained by 
substituting D(z, S\Q) into P^jl - P^]) . When the relationships /x = (3(J m* + h), 7 = 
ppoc*, p = (3 J, and c* = q t are satisfied under the RS assumption, the RS fixed point 
equations (14"T|) - (14"31 are transformed as 

Du tanh (5 (Jom* + h + J^fqu) 



m 



2 V 



J DuDv tanh j3 (Join* + h + J-s/qu) 



x tanh po + Po(a 



Qt 



u + 




(44) 
(45) 

(46) 



In this case, the overlaps are given by 

q = J Du tanh 2 (p + po^/qou) (47) 

q = J .Dm tanh 2 (3 (J m* + h + Jy/qu) , (48) 

and the equations to determine the conjugate po is given by 

m = J Du tanh (p + po^/qou) . (49) 

The conjugates of the energy r and ro are given by p = 2r*/(l — q 2 ) and po = 2ro/ 
(1 — ?q), respectively, and the former relationship indicates the consistency of p = f3J. 
Equations (jUJ, ( H5|) . and fT48|) correspond to the magnetization, randomness energy, 
and overlap in the SK model under the RS assumption, respectively, and hence the RS 
solution of the SK model is recovered as a solution of the fixed point equations, as with 
the original DRT. 



4. Results 



The evolutional dynamics of m, r, and c, given by (!20|) - (|22|) . are solved numerically. 
The continuous evolutional dynamics are treated discretely in the numerical simulation 
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Figure 1. (colour online) Time evolution of (a) m and (b) — r at T/J = 2.0 and 
Jo/ J = 1.0, whose equilibrium state corresponds to a paramagnetic phase. The solid 
and dotted lines represent the dynamics theoretically predicted by the (m, r)-DRT 
and (m,r, c)-DRT, respectively. The dashed line represents the results of the MC 
simulation at N = 4096. 




Figure 2. (colour online) Time evolution of (a) m and (b) — r at T/J = 1.0 and 
Jo / J = 1.5, whose equilibrium state corresponds to a ferromagnetic phase. The lines 
are the same as those in Figure [T] 




2 4 , 6 8 10 2 4 , 6 8 10 



Figure 3. (colour online) Time evolution of (a) m and (b) ~r at T/J = 0.5 and 
Jo/ J — 0.0, whose equilibrium state corresponds to a spin-glass phase. The lines are 
the same as those in Figure [U The short time behaviour of the magnetization is shown 
in the inset of (a). 
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with a unit step of At = 0.01. The changes over time, Am = m(t + At) — m(t), 

Ar = r(t + At) — r(t), and Ac = c(t + At) — c(t), are obtained by calculating the noise 

distribution D(z,S\Q) at each time step t. To do so, the equations of the variables 

(//, // , 7, q, %, qt), f |27|) -( l29|) and ( 132]) - (IM)) . should be solved under the given values of 

the macroscopic quantities fl(t) = {m ,r ,m(t),r(t),c(t)}. The conjugates (/z, ^o,l) 

and overlaps (g, q , q t ) are solved alternately until fl271)- fl29l) and fl32l) - fl34l) are satisfied 

simultaneously; the conjugates are implicitly solved by the Newton method based on 

(I27l)-( l29i) under fixed values of the overlaps, and the overlaps are solved recursively 

based on (I32]) - (I34")) under fixed conjugates. The time dependence of the conjugates and 

overlap parameters are shown in Appendix A , which correspond to the time evolution 

of m, r and c shown in this section. 

The results of the dynamical equations are compared with a Monte Carlo (MC) 

simulation for A" = 4096. The MC simulation settings are as follows. The initial spin 

configuration is randomly generated according to the distribution, 

1 + m 1 ~ra 

p s{Si) = — - — 8 Si ,i + — - — 5 Si ,-i, (50) 

and the interaction J is generated according to the distribution (|2j). The MC dynamics 
are averaged over 1000 realizations of J and 100 initial configurations of S for each J. 
A set of randomly generated spin configurations and interactions generally provides 
r(S) ~ OlN' 1 ), which is sufficiently small to be regarded as zero at iV = 4096. 
Therefore, the average over the initial configuration in the MC simulation corresponds 
to uniform sampling from the subshell characterized by mo and ro = 0, and the overlap 
qo at time t w = is given by m\. In the MC simulation, a randomly chosen fc-th spin is 
flipped with probability Wk given by PJ every 1/N step to match the time scale to that 
of the DRT. 

To clarify the effects of involving the dynamical equation of c into the DRT, the 
dynamics of m and r are also compared with that of the original DRT [8], in which the 
subshell is characterized by m and r only. The original DRT and that including the 
autocorrelation function are denoted by (m, r)-DRT and (to, r, c)-DRT, respectively. 
The results in three parameter regions corresponding to paramagnetic, ferromagnetic, 
and spin-glass phases are shown here. 



4-1. Time evolution of the magnetization and randomness energy 

The time evolution of m and — r are shown in Figures [lj El and [3] for J /J= 1.0 and 
T/J = 2.0 (paramagnetic phase), Jo/ J = 1.5 and T/J = 1.0 (ferromagnetic phase), 
and Jo/ J = 0.0 and T/J = 0.5 (spin-glass phase), respectively. The initial condition at 
t w — is set to m = 0.5 and r = 0. 

In the paramagnetic phase (Figure d]), the time evolution of m and r in the 
MC simulation is well described by the (m, r, c)-DRT proposed here as well as the 
original (to, r)-DRT. The time evolution does not change significantly by including the 
autocorrelation function into the DRT formula in the paramagnetic phase. 
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Figure 4. (colour online) Time evolution of the autocorrelation function c(t w + 1, t w ) 
at Jo/ J = 1.0 and T/J = 2.0. The dashed lines are the results of the MC simulation. 

In the ferromagnetic phase (Figure [2]), the time evolution of m predicted by the 
(m, r, c)-DRT is improved from that of the original (m, r)-DRT. A discrepancy between 
the MC simulation and the DRT results appears in early time steps, but it is confirmed 
that the fixed point of the DRT is coincident with the MC results as t — > oo. The time 
evolution of r described by the (m, r, c)-DRT is almost coincident with that described 
by the original (m, r)-DRT. At sufficiently large t, the DRT predicts the correct value 
of r, but the decay to the fixed point at around 2 < t < 8 is faster than the MC result. 

Figure [3] shows the time evolution of m and r in the spin-glass phase. The short 
time behaviour of m in the (m, r, c)-DRT is slightly improved from that in the (m, r)- 
DRT, as shown in the inset of Figure [3] (a). The fixed point of m in both DRTs is the 
same, i.e. the RS solution of the SK model. The behaviour of r in both DRTs is similar, 
but the long time behaviour in the MC simulation is not correctly predicted by either 
DRT. 

To summarize, by taking into account the time evolution of the autocorrelation 
function, the short time behaviour of the magnetization is improved in the ferromagnetic 
and spin-glass phases, but there is little effect on the randomness energy in any phase. 

4-2. Time evolution of the autocorrelation function 

The time evolution of the autocorrelation function in the (m, r, c)-DRT at a waiting 
time t w is compared to the MC simulation. The autocorrelation c(t w + 1, t w ) for t w > 
is calculated according to ( |T4|) after resetting m = m(t w ), r = r(t w ), and go = q{t w ) 
where m(t w ) and r(t w ) are the magnetization and randomness energy at time t w obtained 
according to the dynamics (112jl and ( TL31) with certain initial conditions mo and tq. 
The overlap at time t w , q(t w ), is determined by ( 152"]) for given values of m(t w ) and 
r(t w ). Therefore, the non-equilibrium state at time t w is assumed to be specified by 
three macroscopic quantities: m(t w ), r(t w ), and q{t w ). The time invariance of the 
autocorrelation function at sufficiently large t w is confirmed over the whole parameter 
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Figure 5. (colour online) Time evolution of the autocorrelation function c(t w + 1, t w ) 
at Jo/ J = 1.5 and T/J = 1.0. The dashed lines are the results of the MC simulation. 

region as c(t) = lim^^.^ c(t w + t,t w ). When one sets r = 0, a solution of the fixed 
point equations is given by p = 7 = 0, as discussed in section 13. 3[ and in this case 
the fixed point of the autocorrelation corresponds to c* = q t = m x m*. Therefore, as 
t — > 00, if this fixed-point solution is stable, the autocorrelation decreases to zero in the 
paramagnetic and spin-glass phases irrespective of the value of mo- On the other hand, 
it converges to a finite value given by m Q m* in the ferromagnetic phase. 

Figure H] shows the time evolution of c(t w + t,t w ) at Jo/ J = 1-0 and Tj J = 2.0, in 
which the equilibrium state corresponds to the paramagnetic phase. Both of the results 
for t w = and t w = 10 show good agreement with the results of the MC simulation. 
The autocorrelation decreases to zero as time increases for any t w . In the paramagnetic 
phase, the solution c* = q t = is a stable solution of the fixed point equations. 

In Figure El the time evolution of c(t w + t, t w ) for t w = and t w = 10 are 
shown for the ferromagnetic phase (Jo/ J = 1.5 and T/J = 1.0). The dynamics of 
the autocorrelation function are well described by the (m, r, c)-DRT for both t w values. 
As the waiting time t w increases, the relaxation of the autocorrelation appears to be 
slow. The converged value of c(t w + t,t w ) corresponds to m{t w ) x m(t w + t) for both 
t w times. In the ferromagnetic phase, c* = qt = mom* is a stable solution of the fixed 
point equations where m* ~ 0.56 in this parameter region. 

In Figure EJ the time evolution of the autocorrelation function in the spin-glass 
phase, Jo/ J = and T/J = 0.5, for t w = and t w = 10 are shown. At t w = 0, the 
short time behaviour of the autocorrelation c is well described by the (m, r, c)-DRT, 
but it deviates from the results of the MC simulation as t increases. The fixed point 
of the autocorrelation function for t w = in the spin-glass phase is c* = 0, as in 
the paramagnetic phase. As t w increases, the discrepancy between the autocorrelation 
function in the (m, r, c)-DRT and that in the MC simulation widens even in short time 
steps. In the spin-glass phase, the initial conditions at t w , mo = m(t w ), r = r(t w ), and 
Qo — q(tw) predicted by the (m, r, c)-DRT, are not coincident with the actual values as 
t w increases, and hence the dynamics after t w in the (m, r, c)-DRT do not agree with 
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(colour online) Time evolution of the autocorrelation function c(t w + 1, t w ) 
and T I J = 0.5. The dashed lines are the results of the MC simulation. 



those in the MC simulation. Generally, DRT predicts a faster relaxation of m and r 
than the actual relaxation as shown in Figure [3j which indicates that DRT describes 
the dynamics at time steps later than the actual time. Therefore, the relaxation of the 
dynamical equation in DRT at t w > is slower than that in the MC simulation. By 
setting the initial conditions of m(t w ) and r(t w ) in the (m, r, c)-DRT to those of the 
actual values in the MC simulation at t w , the short time behaviour of c(t w + t,t w ) is 
improved, although there is still a discrepancy between the results. By comparing the 
results in the paramagnetic and spin-glass phases, it is found that the (m, r, c)-DRT can 
describe the slower relaxation of the autocorrelation function in the spin-glass phase 
than that in paramagnetic phase. 

To characterize the relaxation of the autocorrelation function, it is convenient to 
introduce the integral relaxation time defined by 



The tendency of r int (t w ) as t — > oo is numerically checked by changing the integration 
range t. Figure [7] shows the temperature dependence of Ti nt (t w ) at Jo/ J = and h = 
for t w = 100 with integration ranges t = t w , t = 5t w , and t = 10t w . The spin-glass 
transition temperature in this parameter region is T c j J = 1, as indicated by the vertical 
line in Figured The integral relaxation time in the paramagnetic phase (T / J > 1) does 
not depend on the integration range t, but it increases in the spin-glass phase (T/ J < 1) 
as t increases. This means that T[ nt diverges as t — > oo in the spin-glass phase. Around 
the spin-glass transition temperature, T[ nt exhibits a power law behaviour with respect 
to T/J — T c /J as shown in the inset of Figure [3 In the DRT formulation, T mt goes 
to infinity at T c for T int (t w ) = 1.4(T/ J — T c /J) -06 at sufficiently large t w and t, which 
indicates that the (m, r, c)-DRT can detect the critical slowing down. However, the 
exponent of 0.6 is smaller than the expected value of 2.0 [20]. The proposed (m,r, c)- 
DRT qualitatively describes slow relaxation in the spin-glass phase and the critical 
behaviour around the spin-glass transition temperature, but further development is 
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Figure 7. (colour online) The temperature dependence of the integral relaxation time 
r int at Jq/J = and h — for t w = 100. The vertical line at T c /J = 1.0 indicates the 
spin-glass transition temperature. The inset shows the dependence of Tint on (T — T c )/ 
J with a fitting function denoted by the dashed line. 

necessary for a more accurate description. 
5. Summary and discussion 

In this paper, I have examined the relaxation dynamics in the SK model using DRT and 
including the autocorrelation function. The joint probability of the microscopic states 
at times t w and t w + t was introduced, whose time evolution is described by Glauber 
dynamics. The microscopic states at both times are characterized by macroscopic 
quantities; the magnetization and energy contributed by randomness at time t w + t 
denoted by m and r, those at time t w denoted by mo and ro, as well as the autocorrelation 
function c. Following the DRT procedures, closed dynamical equations were derived 
based on the self-averaging and subshell equipartitioning assumptions. The dynamical 
equations are governed by three overlap parameters q, qo, and q t and the conjugates 
/i,/i , and 7 under the RS assumption. 

The time evolution of m in the (m, r, c)-DRT is improved in the ferromagnetic and 
spin-glass phases compared with (m, r)-DRT but that of r does not change significantly 
in any phase. In the paramagnetic and ferromagnetic phases, the time evolution of 
the autocorrelation function c(t w + 1, t w ) was well described by the proposed framework 
even at a finite t w . In the spin-glass phase, the short time behaviour of c(t w + t,t w ) at 
t w = was well described by the (m, r, c)-DRT, but the theoretical prediction deviated 
from the MC simulation as the waiting time t w and the difference between the two 
times t increased. The temperature dependence of c(t w + t, t w ) was characterized by the 
integral relaxation time r int , and it was found that slow relaxation in the spin-glass phase 
and the critical slowing down are qualitatively described by the current (m, r, c)-DRT 
formulation. 

To improve the results in the spin-glass phase, the validity of the assumptions, RS 
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and subshell equipartitioning with some macroscopic quantities, should be discussed. 
A promising way for the improvement is to introduce a higher-order macroscopic 
quantity to characterize the subshell as in [18] and RSB (TJ |2] . Another straightforward 
development is to take into account multitime correlations. One can obtain the 
dynamical equations of the multitime correlation function by starting with a joint 
probability distribution of the microscopic states at several time steps. In this case, 
the analytical method that is employed to derive the dynamical equations corresponds 
to the replica method for a multi-replicated system consisting of spin systems at each 
time step, which is called real-replica analysis and employed to understand the RSB 
picture in the equilibrium state [2T] . 

DRT has been developed here to include the dynamical equation of the 
autocorrelation function defined by the microscopic states at two specific times. From 
this study, the ability and prospects of DRT for describing the dynamics of random 
systems has been shown. A modification of the DRT proposed in this paper will provide 
an analytical description of the experimentally and numerically employed protocols for 
glassy systems. For example, it is expected that the fluctuation-dissipation relation 
can be checked in the present formulation by comparing the time evolutions of the 
susceptibility and autocorrelation function [22]. The current theory is applicable to a p- 
body interaction system; in this case the relationship between bifurcation in dynamical 
equations and static phase transition associated with RSB should be discussed [23J. 
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Appendix A. Overlap and conjugate parameters 

To obtain the time evolution of the macroscopic variables, the overlaps qo, q, q t and 
conjugates ^0,^,7 should be solved at each time step. Figure lATI and lA"2l show the 
time dependence of the conjugates and overlap parameters, respectively, at t w = 0. 
The initial condition is given by m = 0.5 and r = 0, and the corresponding time 
evolution of the macroscopic quantities are shown in section HJ As shown in Figure IA11 
the conjugate 7 decreases to with increasing time. As t — > 00, the conjugates fi and \x 
converges to the values given by tanh -1 mo and f3(Jom* + h), respectively. The causality 
that q takes a constant value (in this case given by ml) is achieved at all time steps by 
appropriately controlling the conjugates as shown in Figure IA2I The overlap between 
times t w and t w + t, q t , decreases to zero in the paramagnetic and spin-glass phases and 
converges to m x m(t) in the ferromagnetic phase. 
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Figure Al. (colour online) Time dependence of the conjugate parameters at (a) T j 
J = 2.0 and J / J = 1.0, (b) T/J = 1.0 and J /J = 1.5, and (c) T/J = 0.5 and 
Jo/ J = 0.0 for t w = 0. 




Figure A2. (colour online) Time dependence of the overlap parameters at (a) T / 
J = 2.0 and J / J = 1.0, (b) T/J = 1.0 and Jo/ J = 1.5, and (c) T/J = 0.5 and 
Jo/ J = 0.0 for t w = 0. 
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